home *** CD-ROM | disk | FTP | other *** search
/ TOS Silver 2000 / TOS Silver 2000.iso / Anwendun / Pov / POVTOOLS / SOURCE / TRIANGLE.C < prev    next >
Encoding:
C/C++ Source or Header  |  1995-11-25  |  19.4 KB  |  679 lines
  1. /****************************************************************************
  2. *                triangle.c
  3. *
  4. *  This module implements primitives for triangles and smooth triangles.
  5. *
  6. *  from Persistence of Vision Raytracer
  7. *  Copyright 1993 Persistence of Vision Team
  8. *---------------------------------------------------------------------------
  9. *  NOTICE: This source code file is provided so that users may experiment
  10. *  with enhancements to POV-Ray and to port the software to platforms other 
  11. *  than those supported by the POV-Ray Team.  There are strict rules under
  12. *  which you are permitted to use this file.  The rules are in the file
  13. *  named POVLEGAL.DOC which should be distributed with this file. If 
  14. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  15. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  16. *  Forum.  The latest version of POV-Ray may be found there as well.
  17. *
  18. * This program is based on the popular DKB raytracer version 2.12.
  19. * DKBTrace was originally written by David K. Buck.
  20. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  21. *
  22. *****************************************************************************/
  23.  
  24. #include "frame.h"
  25. #include "vector.h"
  26. #include "povproto.h"
  27.  
  28. METHODS Triangle_Methods = 
  29.   {
  30.   All_Triangle_Intersections,
  31.   Inside_Triangle, Triangle_Normal,
  32.   Copy_Triangle,
  33.   Translate_Triangle, Rotate_Triangle,
  34.   Scale_Triangle, Transform_Triangle, Invert_Triangle, Destroy_Triangle
  35.   };
  36.  
  37. METHODS Smooth_Triangle_Methods = 
  38.   {
  39.   All_Triangle_Intersections,
  40.   Inside_Triangle, Smooth_Triangle_Normal,
  41.   Copy_Smooth_Triangle,
  42.   Translate_Smooth_Triangle, Rotate_Smooth_Triangle,
  43.   Scale_Smooth_Triangle, Transform_Smooth_Triangle, 
  44.   Invert_Smooth_Triangle, Destroy_Triangle
  45.   };
  46.  
  47. extern RAY *CM_Ray;
  48. extern long Ray_Triangle_Tests, Ray_Triangle_Tests_Succeeded;
  49.  
  50. #define max3(x,y,z) ((x>y)?((x>z)?1:3):((y>z)?2:3))
  51.  
  52. #define MAX3(x,y,z) (((x)>(y))?(((x)>(z))?(x):(z)):(((y)>(z))?(y):(z)))
  53. #define MIN3(x,y,z) (((x)<(y))?(((x)<(z))?(x):(z)):(((y)<(z))?(y):(z)))
  54.  
  55. void Find_Triangle_Dominant_Axis(Triangle)
  56. TRIANGLE *Triangle;
  57.   {
  58.   DBL x, y, z;
  59.  
  60.   x = fabs(Triangle->Normal_Vector.x);
  61.   y = fabs (Triangle->Normal_Vector.y);
  62.   z = fabs (Triangle->Normal_Vector.z);
  63.   switch (max3(x, y, z)) 
  64.   {
  65.   case 1: 
  66.     Triangle->Dominant_Axis = X_AXIS;
  67.     break;
  68.   case 2: 
  69.     Triangle->Dominant_Axis = Y_AXIS;
  70.     break;
  71.   case 3: 
  72.     Triangle->Dominant_Axis = Z_AXIS;
  73.     break;
  74.   }
  75.   }
  76.  
  77. void Compute_Smooth_Triangle (Triangle)
  78. SMOOTH_TRIANGLE *Triangle;
  79.   {
  80.   VECTOR P3MinusP2, VTemp1, VTemp2;
  81.   DBL x, y, z, uDenominator, Proj;
  82.  
  83.   VSub (P3MinusP2, Triangle->P3, Triangle->P2);
  84.   x = fabs (P3MinusP2.x);
  85.   y = fabs (P3MinusP2.y);
  86.   z = fabs (P3MinusP2.z);
  87.  
  88.   switch (max3 (x, y, z)) 
  89.   {
  90.   case 1:  
  91.     Triangle->vAxis = X_AXIS;
  92.     Triangle->BaseDelta = P3MinusP2.x;
  93.     break;
  94.  
  95.   case 2:  
  96.     Triangle->vAxis = Y_AXIS;
  97.     Triangle->BaseDelta = P3MinusP2.y;
  98.     break;
  99.  
  100.   case 3:  
  101.     Triangle->vAxis = Z_AXIS;
  102.     Triangle->BaseDelta = P3MinusP2.z;
  103.     break;
  104.   }   
  105.  
  106.   VSub (VTemp1, Triangle->P2, Triangle->P3);
  107.   VNormalize (VTemp1, VTemp1);
  108.   VSub (VTemp2, Triangle->P1, Triangle->P3);
  109.   VDot (Proj, VTemp2, VTemp1);
  110.   VScaleEq (VTemp1, Proj);
  111.   VSub (Triangle->Perp, VTemp1, VTemp2);
  112.   VNormalize (Triangle->Perp, Triangle->Perp);
  113.   VDot (uDenominator, VTemp2, Triangle->Perp);
  114.   VInverseScaleEq (Triangle->Perp, -uDenominator);
  115.   }
  116.  
  117. int Compute_Triangle (Triangle,Smooth)
  118. TRIANGLE *Triangle;
  119. int Smooth;
  120.   {
  121.   VECTOR V1, V2, Temp;
  122.   DBL Length, T1, T2, T3;
  123.  
  124.   VSub (V1, Triangle->P1, Triangle->P2);
  125.   VSub (V2, Triangle->P3, Triangle->P2);
  126.   VCross (Triangle->Normal_Vector, V1, V2);
  127.   VLength (Length, Triangle->Normal_Vector);
  128.   /* Set up a flag so we can ignore degenerate triangles */
  129.   if (Length < 1.0e-9)
  130.     {
  131.     Triangle->Degenerate_Flag = TRUE;
  132.     return (0);
  133.     }
  134.  
  135.   /* Normalize the normal vector. */
  136.   VInverseScaleEq (Triangle->Normal_Vector, Length);
  137.  
  138.   VDot (Triangle->Distance, Triangle->Normal_Vector, Triangle->P1);
  139.   Triangle->Distance *= -1.0;
  140.   Find_Triangle_Dominant_Axis(Triangle);
  141.  
  142.   switch (Triangle->Dominant_Axis) 
  143.   {
  144.   case X_AXIS:
  145.     if ((Triangle->P2.y - Triangle->P3.y)*(Triangle->P2.z - Triangle->P1.z) <
  146.       (Triangle->P2.z - Triangle->P3.z)*(Triangle->P2.y - Triangle->P1.y)) 
  147.       {
  148.       Temp = Triangle->P2;
  149.       Triangle->P2 = Triangle->P1;
  150.       Triangle->P1 = Temp;
  151.       if (Smooth) 
  152.         {
  153.         Temp = ((SMOOTH_TRIANGLE *) Triangle)->N2;
  154.         ((SMOOTH_TRIANGLE *) Triangle)->N2 = ((SMOOTH_TRIANGLE *) Triangle)->N1;
  155.         ((SMOOTH_TRIANGLE *) Triangle)->N1 = Temp;
  156.         }
  157.       }
  158.     break;
  159.  
  160.   case Y_AXIS:
  161.     if ((Triangle->P2.x - Triangle->P3.x)*(Triangle->P2.z - Triangle->P1.z) <
  162.       (Triangle->P2.z - Triangle->P3.z)*(Triangle->P2.x - Triangle->P1.x)) 
  163.       {
  164.       Temp = Triangle->P2;
  165.       Triangle->P2 = Triangle->P1;
  166.       Triangle->P1 = Temp;
  167.       if (Smooth) 
  168.         {
  169.         Temp = ((SMOOTH_TRIANGLE *) Triangle)->N2;
  170.         ((SMOOTH_TRIANGLE *) Triangle)->N2 = ((SMOOTH_TRIANGLE *) Triangle)->N1;
  171.         ((SMOOTH_TRIANGLE *) Triangle)->N1 = Temp;
  172.         }
  173.       }
  174.     break;
  175.  
  176.   case Z_AXIS:
  177.     if ((Triangle->P2.x - Triangle->P3.x)*(Triangle->P2.y - Triangle->P1.y) <
  178.       (Triangle->P2.y - Triangle->P3.y)*(Triangle->P2.x - Triangle->P1.x)) 
  179.       {
  180.       Temp = Triangle->P2;
  181.       Triangle->P2 = Triangle->P1;
  182.       Triangle->P1 = Temp;
  183.       if (Smooth) 
  184.         {
  185.         Temp = ((SMOOTH_TRIANGLE *) Triangle)->N2;
  186.         ((SMOOTH_TRIANGLE *) Triangle)->N2 = ((SMOOTH_TRIANGLE *) Triangle)->N1;
  187.         ((SMOOTH_TRIANGLE *) Triangle)->N1 = Temp;
  188.         }
  189.       }
  190.     break;
  191.   }
  192.  
  193.   if (Smooth)
  194.     Compute_Smooth_Triangle((SMOOTH_TRIANGLE *) Triangle);
  195.  
  196.   /* Build the bounding information from the vertices */
  197.   /* Use temps so macro not too big. */
  198.   T1 = MIN3(Triangle->P1.x, Triangle->P2.x, Triangle->P3.x);
  199.   T2 = MIN3(Triangle->P1.y, Triangle->P2.y, Triangle->P3.y);
  200.   T3 = MIN3(Triangle->P1.z, Triangle->P2.z, Triangle->P3.z);
  201.   Make_Vector(&Triangle->Bounds.Lower_Left,T1,T2,T3);
  202.  
  203.   T1 = MAX3(Triangle->P1.x, Triangle->P2.x, Triangle->P3.x);
  204.   T2 = MAX3(Triangle->P1.y, Triangle->P2.y, Triangle->P3.y);
  205.   T3 = MAX3(Triangle->P1.z, Triangle->P2.z, Triangle->P3.z);
  206.   Make_Vector(&Triangle->Bounds.Lengths,T1,T2,T3);
  207.  
  208.   VSub(Triangle->Bounds.Lengths, Triangle->Bounds.Lengths,
  209.     Triangle->Bounds.Lower_Left);
  210.   Triangle->Bounds.Lengths.x += EPSILON;
  211.   Triangle->Bounds.Lengths.y += EPSILON;
  212.   Triangle->Bounds.Lengths.z += EPSILON;
  213.   return (1);
  214.   }
  215.  
  216. int All_Triangle_Intersections (Object, Ray, Depth_Stack)
  217. OBJECT *Object;
  218. RAY *Ray;
  219. ISTACK *Depth_Stack;
  220.   {
  221.   DBL Depth;
  222.   VECTOR IPoint;
  223.  
  224.   if (Intersect_Triangle (Ray, (TRIANGLE *)Object, &Depth))
  225.     {
  226.     VScale (IPoint, Ray -> Direction, Depth);
  227.     VAddEq (IPoint, Ray -> Initial);
  228.     if (Point_In_Clip(&IPoint,Object->Clip))
  229.       {
  230.       push_entry(Depth,IPoint,Object,Depth_Stack);
  231.       return (TRUE);
  232.       }
  233.     }
  234.   return (FALSE);
  235.   }
  236.  
  237. int Intersect_Triangle (Ray, Triangle, Depth)
  238. RAY *Ray;
  239. TRIANGLE *Triangle;
  240. DBL *Depth;
  241.   {
  242.   DBL NormalDotOrigin, NormalDotDirection;
  243.   DBL s, t;
  244.  
  245.   Ray_Triangle_Tests++;
  246.   if(Triangle->Degenerate_Flag)
  247.     return(FALSE);
  248.  
  249.   if (Ray == CM_Ray) 
  250.     {
  251.     if (!Triangle->CMCached) 
  252.       {
  253.       VDot (Triangle->CMNormDotOrigin, Triangle->Normal_Vector, Ray->Initial);
  254.       Triangle->CMNormDotOrigin += Triangle->Distance;
  255.       Triangle->CMNormDotOrigin *= -1.0;
  256.       Triangle->CMCached = TRUE;
  257.       }
  258.  
  259.     VDot (NormalDotDirection, Triangle->Normal_Vector, Ray->Direction);
  260.     if ((NormalDotDirection < Small_Tolerance) &&
  261.       (NormalDotDirection > -Small_Tolerance))
  262.       return (FALSE);
  263.  
  264.     *Depth = Triangle->CMNormDotOrigin / NormalDotDirection;
  265.     }
  266.   else 
  267.     {
  268.     VDot (NormalDotOrigin, Triangle->Normal_Vector, Ray->Initial);
  269.     NormalDotOrigin += Triangle->Distance;
  270.     NormalDotOrigin *= -1.0;
  271.  
  272.     VDot (NormalDotDirection, Triangle->Normal_Vector, Ray->Direction);
  273.     if ((NormalDotDirection < Small_Tolerance) &&
  274.       (NormalDotDirection > -Small_Tolerance))
  275.       return (FALSE);
  276.  
  277.     *Depth = NormalDotOrigin / NormalDotDirection;
  278.     }
  279.  
  280.   if ((*Depth < Small_Tolerance) || (*Depth > Max_Distance))
  281.     return (FALSE);
  282.  
  283.   switch (Triangle->Dominant_Axis) 
  284.   {
  285.   case X_AXIS:
  286.     s = Ray->Initial.y + *Depth * Ray->Direction.y;
  287.     t = Ray->Initial.z + *Depth * Ray->Direction.z;
  288.  
  289.     if ((Triangle->P2.y - s)*(Triangle->P2.z - Triangle->P1.z) <
  290.       (Triangle->P2.z - t)*(Triangle->P2.y - Triangle->P1.y))
  291.       return (FALSE);
  292.  
  293.     if ((Triangle->P3.y - s)*(Triangle->P3.z - Triangle->P2.z) <
  294.       (Triangle->P3.z - t)*(Triangle->P3.y - Triangle->P2.y))
  295.       return (FALSE);
  296.  
  297.     if ((Triangle->P1.y - s)*(Triangle->P1.z - Triangle->P3.z) <
  298.       (Triangle->P1.z - t)*(Triangle->P1.y - Triangle->P3.y))
  299.       return (FALSE);
  300.  
  301.     Ray_Triangle_Tests_Succeeded++;
  302.     return (TRUE);
  303.  
  304.   case Y_AXIS:
  305.     s = Ray->Initial.x + *Depth * Ray->Direction.x;
  306.     t = Ray->Initial.z + *Depth * Ray->Direction.z;
  307.  
  308.     if ((Triangle->P2.x - s)*(Triangle->P2.z - Triangle->P1.z) <
  309.       (Triangle->P2.z - t)*(Triangle->P2.x - Triangle->P1.x))
  310.       return (FALSE);
  311.  
  312.     if ((Triangle->P3.x - s)*(Triangle->P3.z - Triangle->P2.z) <
  313.       (Triangle->P3.z - t)*(Triangle->P3.x - Triangle->P2.x))
  314.       return (FALSE);
  315.  
  316.     if ((Triangle->P1.x - s)*(Triangle->P1.z - Triangle->P3.z) <
  317.       (Triangle->P1.z - t)*(Triangle->P1.x - Triangle->P3.x))
  318.       return (FALSE);
  319.  
  320.     Ray_Triangle_Tests_Succeeded++;
  321.     return (TRUE);
  322.  
  323.   case Z_AXIS:
  324.     s = Ray->Initial.x + *Depth * Ray->Direction.x;
  325.     t = Ray->Initial.y + *Depth * Ray->Direction.y;
  326.  
  327.     if ((Triangle->P2.x - s)*(Triangle->P2.y - Triangle->P1.y) <
  328.       (Triangle->P2.y - t)*(Triangle->P2.x - Triangle->P1.x))
  329.       return (FALSE);
  330.  
  331.     if ((Triangle->P3.x - s)*(Triangle->P3.y - Triangle->P2.y) <
  332.       (Triangle->P3.y - t)*(Triangle->P3.x - Triangle->P2.x))
  333.       return (FALSE);
  334.  
  335.     if ((Triangle->P1.x - s)*(Triangle->P1.y - Triangle->P3.y) <
  336.       (Triangle->P1.y - t)*(Triangle->P1.x - Triangle->P3.x))
  337.       return (FALSE);
  338.  
  339.     Ray_Triangle_Tests_Succeeded++;
  340.     return (TRUE);
  341.   }
  342.   return (FALSE);
  343.   }
  344.  
  345. int Inside_Triangle (IPoint, Object)
  346. VECTOR *IPoint;
  347. OBJECT *Object;
  348.   {
  349.   return (FALSE);
  350.   }
  351.  
  352. void Triangle_Normal (Result, Object, IPoint)
  353. OBJECT *Object;
  354. VECTOR *Result, *IPoint;
  355.   {
  356.   *Result = ((TRIANGLE *)Object)->Normal_Vector;
  357.   }
  358.  
  359. void *Copy_Triangle (Object)
  360. OBJECT *Object;
  361.   {
  362.   TRIANGLE *New;
  363.  
  364.   New = Create_Triangle ();
  365.   *New = * ((TRIANGLE *)Object);
  366.  
  367.   return (New);
  368.   }
  369.  
  370. void Translate_Triangle (Object, Vector)
  371. OBJECT *Object;
  372. VECTOR *Vector;
  373.   {
  374.   TRIANGLE *Triangle = (TRIANGLE *) Object;
  375.   VECTOR Translation;
  376.  
  377.   VEvaluate (Translation, Triangle->Normal_Vector, *Vector);
  378.   Triangle->Distance -= Translation.x + Translation.y + Translation.z;
  379.   VAddEq (Triangle->P1, *Vector)
  380.     VAddEq (Triangle->P2, *Vector)
  381.       VAddEq (Triangle->P3, *Vector)
  382.  
  383.         /* Recalculate the bounds */
  384.         VAddEq(Object->Bounds.Lower_Left, *Vector);
  385.   }
  386.  
  387. void Rotate_Triangle (Object, Vector)
  388. OBJECT *Object;
  389. VECTOR *Vector;
  390.   {
  391.   TRANSFORM Trans;
  392.  
  393.   Compute_Rotation_Transform (&Trans, Vector);
  394.   Transform_Triangle (Object, &Trans);
  395.   }
  396.  
  397. void Scale_Triangle (Object, Vector)
  398. OBJECT *Object;
  399. VECTOR *Vector;
  400.   {
  401.   TRIANGLE *Triangle = (TRIANGLE *) Object;
  402.   DBL Length,T1,T2,T3;
  403.  
  404.   Triangle->Normal_Vector.x = Triangle->Normal_Vector.x / Vector->x;
  405.   Triangle->Normal_Vector.y = Triangle->Normal_Vector.y / Vector->y;
  406.   Triangle->Normal_Vector.z = Triangle->Normal_Vector.z / Vector->z;
  407.  
  408.   VLength(Length, Triangle->Normal_Vector);
  409.   VInverseScaleEq (Triangle->Normal_Vector, Length);
  410.   Triangle->Distance /= Length;
  411.  
  412.   VEvaluateEq (Triangle->P1, *Vector);
  413.   VEvaluateEq (Triangle->P2, *Vector);
  414.   VEvaluateEq (Triangle->P3, *Vector);
  415.  
  416.   /* Recompute the bounds */
  417.   /* Use temps so macro not too big. */
  418.   T1 = MIN3(Triangle->P1.x, Triangle->P2.x, Triangle->P3.x);
  419.   T2 = MIN3(Triangle->P1.y, Triangle->P2.y, Triangle->P3.y);
  420.   T3 = MIN3(Triangle->P1.z, Triangle->P2.z, Triangle->P3.z);
  421.   Make_Vector(&Triangle->Bounds.Lower_Left,T1,T2,T3);
  422.  
  423.   T1 = MAX3(Triangle->P1.x, Triangle->P2.x, Triangle->P3.x);
  424.   T2 = MAX3(Triangle->P1.y, Triangle->P2.y, Triangle->P3.y);
  425.   T3 = MAX3(Triangle->P1.z, Triangle->P2.z, Triangle->P3.z);
  426.   Make_Vector(&Triangle->Bounds.Lengths,T1,T2,T3);
  427.  
  428.   VSub(Triangle->Bounds.Lengths, Triangle->Bounds.Lengths,
  429.     Triangle->Bounds.Lower_Left);
  430.   Triangle->Bounds.Lengths.x += EPSILON;
  431.   Triangle->Bounds.Lengths.y += EPSILON;
  432.   Triangle->Bounds.Lengths.z += EPSILON;
  433.   }
  434.  
  435. void Transform_Triangle (Object, Trans)
  436. OBJECT *Object;
  437. TRANSFORM *Trans;
  438.   {
  439.   TRIANGLE *Triangle = (TRIANGLE *) Object;
  440.  
  441.   MTransPoint (&Triangle->Normal_Vector,
  442.     &Triangle->Normal_Vector, Trans);
  443.   MTransPoint (&Triangle->P1, &Triangle->P1, Trans);
  444.   MTransPoint (&Triangle->P2, &Triangle->P2, Trans);
  445.   MTransPoint (&Triangle->P3, &Triangle->P3, Trans);
  446.   Compute_Triangle (Triangle,FALSE);
  447.   }
  448.  
  449. TRIANGLE *Create_Triangle()
  450.   {
  451.   TRIANGLE *New;
  452.  
  453.   if ((New = (TRIANGLE *) malloc (sizeof (TRIANGLE))) == NULL)
  454.     MAError ("triangle");
  455.  
  456.   INIT_OBJECT_FIELDS(New,TRIANGLE_OBJECT,&Triangle_Methods)
  457.  
  458.     Make_Vector (&(New -> Normal_Vector), 0.0, 1.0, 0.0);
  459.   New -> Distance = 0.0;
  460.   New -> CMNormDotOrigin = 0.0;
  461.   New -> CMCached = FALSE;
  462.   Make_Vector (&(New -> P1), 0.0, 0.0, 0.0);
  463.   Make_Vector (&(New -> P2), 1.0, 0.0, 0.0);
  464.   Make_Vector (&(New -> P3), 0.0, 1.0, 0.0);
  465.   New -> Degenerate_Flag = FALSE;
  466.  
  467.   /* NOTE: Dominant_Axis is computed when Parse_Triangle calls
  468.    Compute_Triangle.  vAxis is used only for smooth triangles */
  469.  
  470.   return (New);
  471.   }
  472.  
  473. void Invert_Triangle (Object)
  474. OBJECT *Object;
  475.   {
  476.   return;
  477.   }
  478.  
  479. /* Calculate the Phong-interpolated vector within the triangle
  480.    at the given intersection point. The math for this is a bit
  481.    bizarre:
  482.  
  483.     -         P1
  484.     |        /|\ \
  485.     |       / |Perp\
  486.     |      /  V  \   \
  487.     |     /   |    \   \
  488.   u |    /____|_____PI___\
  489.     |   /     |       \    \
  490.     -  P2-----|--------|----P3
  491.               Pbase    PIntersect
  492.         |-------------------|
  493.                        v
  494.  
  495.    Triangle->Perp is a unit vector from P1 to Pbase. We calculate
  496.  
  497.    u = (PI - P1) DOT Perp / ((P3 - P1) DOT Perp).
  498.  
  499.    We then calculate where the line from P1 to PI intersects the line P2 to P3:
  500.    PIntersect = (PI - P1)/u.
  501.  
  502.    We really only need one coordinate of PIntersect.  We then calculate v as:
  503.  
  504.       v = PIntersect.x / (P3.x - P2.x)
  505.  or   v = PIntersect.y / (P3.y - P2.y)
  506.  or   v = PIntersect.z / (P3.z - P2.z)
  507.  
  508.    depending on which calculation will give us the best answers.
  509.  
  510.    Once we have u and v, we can perform the normal interpolation as:
  511.  
  512.      NTemp1 = N1 + u(N2 - N1);
  513.      NTemp2 = N1 + u(N3 - N1);
  514.      Result = normalize (NTemp1 + v(NTemp2 - NTemp1))
  515.  
  516.    As always, any values which are constant for the triangle are cached
  517.    in the triangle.
  518. */
  519.  
  520. void Smooth_Triangle_Normal (Result, Object, IPoint)
  521. OBJECT *Object;
  522. VECTOR *Result, *IPoint;
  523.   {
  524.   SMOOTH_TRIANGLE *Triangle = (SMOOTH_TRIANGLE *) Object;
  525.   VECTOR PIMinusP1, NTemp1, NTemp2;
  526.   DBL u = 0.0, v = 0.0;
  527.  
  528.   VSub (PIMinusP1, *IPoint, Triangle->P1);
  529.   VDot (u, PIMinusP1, Triangle->Perp);
  530.   if (u < 1.0e-9) 
  531.     {
  532.     *Result = Triangle->N1;
  533.     return;
  534.     }
  535.  
  536.   /* BaseDelta contains P3.x-P2.x,  P3.y-P2.y, or P3.z-P2.z depending on the
  537.       value of vAxis. */
  538.  
  539.   switch (Triangle->vAxis) 
  540.   {
  541.   case X_AXIS:  
  542.     v = (PIMinusP1.x/u + Triangle->P1.x - Triangle->P2.x) / Triangle->BaseDelta;
  543.     break;
  544.  
  545.   case Y_AXIS:  
  546.     v = (PIMinusP1.y/u + Triangle->P1.y - Triangle->P2.y) / Triangle->BaseDelta;
  547.     break;
  548.  
  549.   case Z_AXIS:  
  550.     v = (PIMinusP1.z/u + Triangle->P1.z - Triangle->P2.z)/ Triangle->BaseDelta;
  551.     break;
  552.   }
  553.  
  554.   VSub (NTemp1, Triangle->N2, Triangle->N1);
  555.   VScaleEq (NTemp1, u);
  556.   VAddEq (NTemp1, Triangle->N1);
  557.   VSub (NTemp2, Triangle->N3, Triangle->N1);
  558.   VScaleEq (NTemp2, u);
  559.   VAddEq (NTemp2, Triangle->N1);
  560.   VSub (*Result, NTemp2, NTemp1);
  561.   VScaleEq (*Result, v);
  562.   VAddEq (*Result, NTemp1); 
  563.   VNormalize (*Result, *Result);
  564.   }
  565.  
  566. void *Copy_Smooth_Triangle (Object)
  567. OBJECT *Object;
  568.   {
  569.   SMOOTH_TRIANGLE *New;
  570.  
  571.   New = Create_Smooth_Triangle ();
  572.   *New = * ((SMOOTH_TRIANGLE *)Object);
  573.  
  574.   return (New);
  575.   }
  576.  
  577. void Translate_Smooth_Triangle (Object, Vector)
  578. OBJECT *Object;
  579. VECTOR *Vector;
  580.   {
  581.   SMOOTH_TRIANGLE *Triangle = (SMOOTH_TRIANGLE *) Object;
  582.   VECTOR Translation;
  583.  
  584.   VEvaluate (Translation, Triangle->Normal_Vector, *Vector);
  585.   Triangle->Distance -= Translation.x + Translation.y + Translation.z;
  586.   VAddEq (Triangle->P1, *Vector)
  587.     VAddEq (Triangle->P2, *Vector)
  588.       VAddEq (Triangle->P3, *Vector)
  589.         Compute_Triangle ((TRIANGLE *) Triangle,TRUE);
  590.   }
  591.  
  592. void Rotate_Smooth_Triangle (Object, Vector)
  593. OBJECT *Object;
  594. VECTOR *Vector;
  595.   {
  596.   TRANSFORM Trans;
  597.  
  598.   Compute_Rotation_Transform (&Trans, Vector);
  599.   Transform_Smooth_Triangle (Object, &Trans);
  600.   }
  601.  
  602. void Scale_Smooth_Triangle (Object, Vector)
  603. OBJECT *Object;
  604. VECTOR *Vector;
  605.   {
  606.   SMOOTH_TRIANGLE *Triangle = (SMOOTH_TRIANGLE *) Object;
  607.   DBL Length;
  608.  
  609.   Triangle->Normal_Vector.x = Triangle->Normal_Vector.x / Vector->x;
  610.   Triangle->Normal_Vector.y = Triangle->Normal_Vector.y / Vector->y;
  611.   Triangle->Normal_Vector.z = Triangle->Normal_Vector.z / Vector->z;
  612.  
  613.   VLength(Length, Triangle->Normal_Vector);
  614.   VScaleEq (Triangle->Normal_Vector, 1.0 / Length);
  615.   Triangle->Distance /= Length;
  616.  
  617.   VEvaluateEq (Triangle->P1, *Vector);
  618.   VEvaluateEq (Triangle->P2, *Vector);
  619.   VEvaluateEq (Triangle->P3, *Vector);
  620.   Compute_Triangle ((TRIANGLE *) Triangle,TRUE);
  621.   }
  622.  
  623. void Transform_Smooth_Triangle (Object, Trans)
  624. OBJECT *Object;
  625. TRANSFORM *Trans;
  626.   {
  627.   SMOOTH_TRIANGLE *Triangle = (SMOOTH_TRIANGLE *) Object;
  628.  
  629.   MTransPoint (&Triangle->Normal_Vector,
  630.     &Triangle->Normal_Vector, Trans);
  631.   MTransPoint (&Triangle->P1, &Triangle->P1, Trans);
  632.   MTransPoint (&Triangle->P2, &Triangle->P2, Trans);
  633.   MTransPoint (&Triangle->P3, &Triangle->P3, Trans);
  634.   MTransPoint (&Triangle->N1, &Triangle->N1, Trans);
  635.   MTransPoint (&Triangle->N2, &Triangle->N2, Trans);
  636.   MTransPoint (&Triangle->N3, &Triangle->N3, Trans);
  637.   Compute_Triangle ((TRIANGLE *) Triangle,TRUE);
  638.   }
  639.  
  640. void Invert_Smooth_Triangle (Object)
  641. OBJECT *Object;
  642.   {
  643.   return;
  644.   }
  645.  
  646. SMOOTH_TRIANGLE *Create_Smooth_Triangle()
  647.   {
  648.   SMOOTH_TRIANGLE *New;
  649.  
  650.   if ((New = (SMOOTH_TRIANGLE *) malloc (sizeof (SMOOTH_TRIANGLE))) == NULL)
  651.     MAError ("smooth triangle");
  652.  
  653.   INIT_OBJECT_FIELDS(New,SMOOTH_TRIANGLE_OBJECT,&Smooth_Triangle_Methods)
  654.  
  655.     Make_Vector (&(New->Normal_Vector), 0.0, 1.0, 0.0);
  656.   New->Distance = 0.0;
  657.   New -> CMNormDotOrigin = 0.0;
  658.   New -> CMCached = FALSE;
  659.   Make_Vector (&(New -> P1), 0.0, 0.0, 0.0);
  660.   Make_Vector (&(New -> P2), 1.0, 0.0, 0.0);
  661.   Make_Vector (&(New -> P3), 0.0, 1.0, 0.0);
  662.   Make_Vector (&(New -> N1), 0.0, 1.0, 0.0);
  663.   Make_Vector (&(New -> N2), 0.0, 1.0, 0.0);
  664.   Make_Vector (&(New -> N3), 0.0, 1.0, 0.0);
  665.   New -> BaseDelta = 0.0;
  666.   New -> Degenerate_Flag = FALSE;
  667.  
  668.   /* NOTE: Dominant_Axis and vAxis are computed when Parse_Triangle calls
  669.    Compute_Triangle.  */
  670.  
  671.   return (New);
  672.   }
  673.  
  674. void Destroy_Triangle (Object)
  675. OBJECT *Object;
  676.   {
  677.   free (Object);
  678.   }
  679.